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Tree-Structured Grid Model of Line and Polarization Variability 

from Massive Binaries 
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Abstract. We have developed a 3-D Monte Carlo radiative transfer model which computes line and continuum 
polarization variability for a binary system with an optically thick non-axisymmetric envelope. This allows us 
to investigate the complex (phase-locked) line and continuum polarization variability features displayed by many 
massive binaries: W-R+O, O+O, etc. An 8-way tree data structure constructed via a "cell-splitting" method 
allows for high precision with efficient use of computer resources. The model is not restricted to binary systems; 
it can easily be adapted to a system with an arbitrary density distribution and large density gradients. As an 
application to a real system, the phase dependent Stokes parameters (I, Q, U) and the phase dependent He I 
(A5876) profiles of the massive binary system V444 Cyg (WN5+06 III-V) are computed. 
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1. Introduction 

Linear polarization of light is often used to study the geo- 
metrical distribution of gaseous/dusty materials around 
various types of astronomical objects, e.g., close bina- 
ries, Be stars, Wolf-Rayet (W-R) stars, Luminous Blue 
Variables (LBVs), supernovae, Seyfert galaxies, quasars 
(QSOs), and Young Stellar Objects (YSOs). When the 
angular size of an object is too small to be resolved by 
ground-based or space telescopes, spectropolarimetry and 
narrow band polarimetry provide useful information on 
the geometry of circumstellar matter. 

Several analytic models use simplifying assumptions to 
predict the continuum polarization levels arising from elec- 
tron scattering envelopes. For example, Brown & McLean 
|(1977D derived an expression for the degree of polariza- 
tion produced in an optically thin axisymmetric envelope 
with an embedded point source. Brown et aL (1978) ex- 
tended this model, deriving the expressions for the polar- 
ization arising from a binary system embedded in an op- 
tically thin envelope. Adopting this model, Brown et al 
|(1989D considered the effect of the finite size of the light 
source. Unfortunately, these 2-D models lack an ability 
to investigate more complex systems such as close bina- 
ries with colliding stellar winds, molecular clouds with 
fractal structure, and the tails of non-circular accretion 
flows. Further, they can not deal with the realistic in- 
homogeneities present in many scattering envelopes, and 



they are not suitable for optically thick atmospheres with 
multiple scattering. 

In order to investigate such complicated systems, we 
abandon the simplifying assumptions and use a full three 
dimensional model that can handle an optically thick at- 
mosphere. Formulating and solving the radiative transfer 
equation for an arbitrary geometry is difficult even for 
continuum radiation. The presence of additional scatter- 
ing integrals in the expression for the source function and 
the need to solve several coupled transfer equations for 
the Stokes parameters tremendously complicates the prob- 
lem. The flexibility of Monte Carlo methods of radiative 
transfer are particularly useful for such problems despite 
the fact that they may not be the most efficient solution 
method. For example, asymmetries in the radiation field, 
non-spherical geometry, and a complex distribution of the 
clumps can be readily incorporated. The issue of efficiency 
is less severe due to the significant advancement of com- 
puter technology. In the near future, it may be possible 
to perform non-LTE radiative transfer calculations in 3- 
D including thousands of lines from m ulti-atomic sp ecies 
utilizing some of the i deas describ ed by Bernes (1979 ) and 
Li fc McCray (1996[ ). Recently, |Lucy (1999| ) developed 
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an improved Monte Carlo technique for non-LTE radia- 
tive transfer calculations applied to a synthetic supernova 
spectra model. 

Our main reason for developing a 3-D model is to 
investigate the colliding wind (CW) interaction zone. 
Such zones are relatively common among massive bi- 
naries (W-R+O, O+O binaries). |Koenigsberger fc Auei 
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(1985) found, using low resolution IUE spectra, that 3 
out of 6 selected W-R+O binaries showed evidence for 
CWs. [Shore fc Brown (19881 ), |Marchenko et al. (1994| ) and 
Marchenko ct al. (1997) presented detailed studies of CW 
related variability seen in high-resolution UV and optical 
data of V444 Cyg (WN5+06 III-V). |5tevens fc Howarth 



(1999|)]presented observations of the phase dependent He I 
(1.0830-/mi) line profile in six different W-R binary sys- 
tems including V444 Cyg. Using a simple model, they ex- 
plained the observed variability in terms of the wind- wind 
interaction. Further theoretical aspects of the relationship 
bet ween X-ray varia bility and CWs have been developed 
by |Luo et al. (1990|) and [Stevens ct al. (1992|) . More re- 



cently, detailed hydrodynamic models of c olliding winds 
have be use d by [Stevens fc Pollock (1994 ), Gayley et al 



(19970 and |Pittard (1998j ). Si mulations of the collidin g 
winds in V444 Cyg are given in [Pittard fc Stevens (1999| ) 



Examples of 3-D Mont e Carlo radiative trans f er mod- 
els rece ntly developed are Witt fc Gordon (1996 ), Pagani 
(1998|)] [Stevens fc Howarth (1999| ), |Wolf et al. (19991 ) and 
Harries (200C ). Juvela (1997 ) presented their 1-3 dimen- 



sional Monte Carlo models of emission from clumpy molec- 
ular clouds. Wolf et al. (1999| ) developed a self-consistent 
3-D continuum Monte Carlo model in which the dust tem- 
perature and the radiation field are iterated to consis- 



tency. These 3-D models, with the exceptions of Wolf et al 
(1999|)]and Harries (200C| ), use regular "cubic" grids and 



are not suitable for a model with a large density gradient 
or a large dynamic range of density. These would require 
much finer grid sizes (and hence more memory) to re- 
solve the smaller structures. |Wolf et al. (1999] ) and Harries 



are evenly spaced in polar and azimuthal angles, but log- 
arithmically spaced in the radial direction. This type of 
grid is useful for a system that does not deviate much 
from spherical or axisymmetric symmetry, but it would 
not be adequate for a binary system with colliding stellar 
winds or for a binary system with accretion flow from a 
companion 

Starting from the 2-D Monte Carlo model of |Hillicr 



(199l|)] we have developed a 3-D Monte Carlo code to 
calculate the continuum and line polarizations produced 
by the scattering of light in an arbitrary geometry. The 
model can predict the variability features associated with 
the orbital motion of a binary system, such as the polar- 
ization level, flux level, and line profile shapes. The model 
can treat a finite size stellar disk, multiple scattering, ab- 
sorption of continuum photons by a line, and emission 
from multiple light sources (extended or not) . Higher pre- 
cision is achieved with fewer grid points by using a "cell- 
splitting" method whose application to this type of prob- 
lem was discussed by Wolf et al. (1999| ) . We extended the 
idea of "cell-splitting" to a simple 8-way tree data struc- 
ture to construct the model grids. With this data struc- 
ture, we can quickly access the data stored at any grid 
point, leading to a faster, more efficient, numerical code. 
This is essential to our model because of the high dimen- 
sionality and because no assumptions are made regarding 



the symmetry of the model. A logarithmic grid, like those 
commonly used in spherical and axisymmetric codes, is 
not readily implementable. 

In § J2| we describe the details of our model. Various 
accuracy and efficiency tests are given in § |. The model 
is applied to the massive binary V444 Cyg (WN5+06) 
as an example calculation for a real system in § ^. The 
conclusions are given in § ||. 

2. Model 

2.1. Polarization 

A general method of polarization calculation via the 



Monte Carlo method is described in Modali et al. (1972), 
|Sandford (1973] ), |Warren-Smith (1983| ) and [Hillier (199l| ). 
An initial photon beam emitted from a light source is 
usually assumed to be unpolarized, with the Stokes pa- 
rameters given by (/, Q,U, V) = (I , 0,0,0). After each 
scattering, the Stokes parameters change according to the 
Mueller matrix, M(Q, <p) where 9 and <f> are the polar and 
azimuthal scattering angles respectively. 
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In the case of Thomson scattering , M has no 4> depen- 
dency, and it becomes 



(2000j)]used 3-D grids in spherical coordinates. Their grids M 
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where (i — cos 9. 

A coordinate transformation of the Stokes parame- 
ter vector A T = (/, Q, U, V) to a fixed reference frame 
is necessary after each scattering. Chandrasekhar (1960) 
describes how this transformation is done. After several 
scattering events, the photon beam will reach the model 
boundary (cubic in our case), and the A T vector will then 
be proj ected onto th e plane of an observer. The reader is 
refer to Hillier (1991 ) for details. 



2.2. Coordinates 

For the model discussed here, we assume that the density 
distribution is co-rotating with the binary system, and the 
orbit is circular. Two coordinate systems are used in the 
model (see Fig. |l]) . One is the primary star coordinates (x, 
y, z), and the other is the binary coordinates (x', y', z'). 
The binary coordinate system is a rotating coordinate sys- 
tem, and it is chosen such that the stars and the density 
are fixed (or independent of orbital phases) in that coor- 
dinate system. The orbital plane is assigned to be on the 
x-y plane of the primary star coordinates. The primary 
star (Star A) is placed at the origin of the two coordi- 
nate systems, and the secondary star (Star B) is placed 
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z, z 



To observer 




Fig. 1. This illustrates the two coordinate systems used 
in the models. The binary coordinates (x', y' z') rotates 
around the primary star coordinates (x, y, z) through the 
common z or z' axis. A rotation angle (0) and a binary 
separation (a) are functions of orbital phase t. For a circu- 
lar orbit, they are simply, 8 (t) = 2n and a =constant. An 
observer is placed on the y-z plane with the viewing angle 
i. The primary star (Star A) is placed at the origin of the 
two coordinate systems, and the secondary star (Star B) 
is placed on the x' axis. With this configuration, the orbit 
of the binary is confined on x-y plane. 

on the x' axis. The relative position of two coordinate 
systems changes as a function of orbital phase (t). When 
t = 0, the two coordinate systems coincide. For t > 0, the 
primed (binary) coordinate system rotates around the z 
axis of the unprimed (primary star) coordinates making 
the angle between x & x' and y & y' to be (t). For a 
circular orbit the relative angle is simply, (t) = 2nt. In 
general, the separation distance a is a function of t, but 
for a circular orbit a remains constant. By our convention, 
an observer is always on the y-z plane with an inclination 
angle i measured from z or z' axis. 

In the primary star (unprimed) coordinates, the loca- 
tions of the stars, as a function of i, become 



x A (t) 

y A (t) 

z A {t) 

and 

x B {t) 

y B (t) 



(1) 



a cos {2n (t + to)} 
a sin {2tt (t + to)} 




(2) 



where 2irt is the relative angle between the two coor- 
dinates and to is an initial phase offset. In our model 
to = —1/4 is used so that when t = 0, the secondary 
star (Star B) is behind the primary star (Star A) from an 
observer's perspective. Eqs. |l| andg can be easily modified 



to an orbit with non-zero eccentricity for a general binary 
orbit. 



2.3. Grid Selection 

When the gradient in the opacity field is very large, a 
logarithmic scale in the radial direction can be used to in- 
crease computational accuracy of the optical depth if the 
system is spherically or axially symmetric. For the case 
of an arbitrary geometry, there is no simple way to con- 
struct a logarithmic scale. However, an efficient cubic grid 
(consisting of cubic cells) can be constructed by subdi- 
viding a cube into smaller cubes where the value of the 
opacity /cmissivity is larger than a given threshold value. 



Wolf et al. (1999) introduced the basic idea for this type of 



grid method to the radiative transfer problem. We do not 
know the exact algorithm they used, but we have utilized 
the 8-way tree data structure to construct the grids and 
to store the values of the opacities. In addition to being 
conceptually simple, the advantages of using a tree data 
structure include increases in optical depth calculation ac- 
curacy, computer memory efficiency, and data access time. 
In the following subsections, we describe the algorithm for 
constructing the tree-structured grid and illustrate its ef- 
ficiency. 

2.4. Tree Data Structure Construction Algorithm 

The data structure used here is similar to a "binary tree" 
in which one node splits into two children nodes, but ours 
splits into eight. The steps for the tree construction are 
summarized below. Fig. || gives the flow chart of the steps. 

1. Start from a cubic box {root node) which holds the 
whole density structure. 

2. Pick (100-1000) random locations, r = (x,y,z), in the 
cell, and compute 



E 



n 



(3) 



where r\i is the emissivity at i-th random location, n is 
the total number of the random locations, d is the size 
of the cell and p is the index of scalability. 

3. If E is greater than E Q which is a user specified pa- 
rameter, divide the box into 8 equal-size cubes (child 
nodes). If E is smaller than E Q , do nothing. 

4. For each cell just created, repeat Steps 2 and 3 until all 
cells have E < E Q . — The lowest level cells are often 
called "leaf(s)" (leaves). 

The index p in Eq. || controls how fast the cell should be 
split with 77, and normally p = 1 is used. When p = 1, 
E represents the average "emission measure" of the cell 
since it simply becomes the product of the average emis- 
sivity and the volume of the cell. The emissivity at a given 
point can be evaluated from one or a combination of the 
following: l.an analytical formula, 2. an interpolation of 
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Fig. 2. This diagram shows how a cell in a model space 
will be split into 8 subcells. 1. Consider a cubic cell with 
length "d." 2. Choose n random locations in the cell, and 
evaluate the emissivity value (rji) at each location. Average 
these emissivities to obtain E. 3. If E is greater than a 
threshold value E a , then 4. Divide the cell into 8 equal 
size cubic cells. The program will repeat this procedure 
recursively until all the cells satisfy the condition: E < E Q . 
Note that the figure shows the cell split into only 4 cells 
because the 3-D cell is projected onto a 2-D plane for 
clarity. 



the output from a radiative transfer model and 3. an in- 
terpolation of the output from a hydrodynamical model. 
Other more complex algorithms could be used to optimize 
the computation for a specific problem. 

Once the data structure is built, the access/search for 
a given cell can be done recursively, resulting a faster code. 
Examples of the resulting grids are shown in Fig. |§|. The 
top diagram in Fig. || shows the grids constructed for two 
stars with different sizes having spherically distributed gas 
around each star. The figure clearly shows that the grid 
size becomes naturally smaller where the density is higher. 
The grids used for this density is in 3-D, but the diagram 
shows only the cross section of the y-z plane for clarity. In 
the bottom diagram of Fig. ||, the grids are overlaid with a 
pseudo-spiral galaxy density distribution (3-D but drawn 
in 2-D). The density of the disk is inversely proportional 
to the distance from the center of the galaxy except for 
the bulge region where the density is constant. The grids 
nicely follow the spiral arms, and become smaller as they 
approach to the center. 



2.5. Searching for a Cell 

An example of how to find the cell which contains a point 
(p) in the tree data structure is illustrated in Fig. ||. The 
steps shown here are for the 2-D case for clarity, and the 





Fig. 3. Examples of the tree structured grids. Top: a bi- 
nary system with different stellar radii surrounded by a 
spherical atmosphere around each star. Bottom: a two- 
arm pseudo spiral galaxy. The density decreases as oc 1 /r 2 
for r > Rb and it is constant for a r < Rb where Rb is the 
bulge radius. Note that the real density distributions are 
3-D, but here they are depicted in 2-D for simplicity. In 
both cases, the grid size naturally decreases for the high 
density regions. 



same method can be applied to the 3-D case. The figure 
shows the following steps: 

1. Starting from the root cell (the outer most cell), out- 
lined with a thick solid line, check whether the point 
is in the left half or the right half of the cell. 

2. The point is in the right half of the cell; hence, dismiss 
all the subcells in the left half. 

3. Further, check whether the point is in the upper half 
or the lower half of the cell outlined with a thick solid 
line. The point is in the upper half of the cell; hence, 
dismiss all the subcells in the lower half. 
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Fig. 4. An illustration of how the cell which contains a 
point (p) in the model space can be found in the tree data 
structure. The figure demonstrates the searching steps for 
the 2-D case for clarity, but the same method can be ap- 
plied to the 3-D case. The thick solid line indicates the 
cell that is focused in each step. The thick dotted lines 
represent the border used to check which side of the cell 
point p occupies. The corresponding instructions for each 
step are summarized in 
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4. Further, check whether the point is in the left half or 
the right half of the cell outlined with a thick solid line. 

5. The point is in the right half of the cell; hence, dismiss 
all the subcells in the left half. Further, check whether 
the point is in the upper or the lower of the cell outlined 
with a thick solid line. 

6. The point is in the lower half of the cell. Since this is 
the lowest level (leaf node) of the tree, the search stops 
here. The leaf cell containing point p is found. 

In our algorithm, whether a point belongs to the left 
half or the right half of the cell and whether a point be- 
longs to the upper half or the lower half of the cell are 
checked simultaneously. Therefore, steps 1 and 2 are con- 
sidered to be just one step, as are steps 3 and 4. In other 
words, the total number of steps required to find the cell 
including point p is considered to be "two" in this exam- 
ple. 



2.6. Performance Check 

To demonstrate the efficiency of the tree structure |] the 
optical depth and the volume integral of density were 
computed with a (single-size) regular cubic grid and a 
tree structure grid separately. The functional forms of 
the opacity and the density used are x ( r ) — Xo/j 4 an d 
p (r) — po/r 2 respectively. The top graph in Fig. shows 
the percentage error of the optical depth calculations with 
the tree structure grid and that with one-size cubic grids 
as a function of "total" number of grid points. Similarly, 
the bottom graph in Fig. |B| shows the percentage error 
of volume integrals as a function of total number of grid 
points. Note that no special integration weights were used 
in both cases. 

The figure shows the error decreases very quickly as 
the total number of grid/node points increases for the tree 
structure grids. As a consequence, much smaller number of 
points is required for the tree method to achieve the same 
order of accuracy, and thus the method requires much less 
computer memory. This enables us to do a realistic 3-D 
simulation even on a regular PC with a reasonable amount 
of RAM. 

The number of steps in a grid searching routine de- 
pends on the complexities of the data structure, but in 
the case of the grids constructed with the inverse-square 
density function with the total grid points of ~ 10 6 , the 
average number of searching steps is about 5. Although 
not implemented to the current model, the searching could 
be made more efficient by searching upward from the cur- 
rent position of a photon rather than searching downward 
from the root node each time since the cells along a photon 
trajectory are correlated. 

2. 7. Optical Depth 

Suppose a photon at r Q = (x y ,z ) is moving in the 
direction of N as shown in Fig. ^. This photon will most 
likely intersect with many cells (or grids) before it reaches 
the outer boundary. Calculation of the optical depth for a 
photon which is at r D and moving in direction N is done 
by replacing the integral along the ray with the sum of the 
product of the cell opacity {xi) with the line length (8si) 
of the photon beam in the cell, yielding 



X(s) ds 



(4) 



where Xi is the sum of thermal and electron scattering 
opacities for a continuum photon: 



Xi 



Xi ' Xi 



and m is the number of cells which intersect with the pho- 
ton. In order to perform the summation in Eq. m we first 



1 For this simple density distribution, the cylindrical coor- 
dinates (r, 4>, z) with logarithmic spacing in r direction would 
be very efficient, but such a grid scheme is too restricted for 
an arbitrary density distribution. 
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Optical Depth Integration 




number of grid points 



Fig. 5. Top: Percentage errors of optical depth integrals 
using a regular grid (squares) and using a tree structured 
grid (circles) v.s. total number of grid points. The opacity 
is set to be \ — Xo/r 4 , and the optical depth is calcu- 
lated from the surface of a star (r = li?*) to the outer 
boundary of the model (r = 100i?»). Bottom: Percentage 
errors of volume integrals (J2aii ceils Pi^Vi) using a regular 
grid (squares) and using a tree structured grid (circles) 
v.s. total number of grid points. The density is assumed 
to be p = p /f 2 , the integration limits are r min = li?* and 
'"max = 100i?* where i?» is the radius of a star. 



need to find with which cells the photon will intersect. 
Then, the opacity values stored in these cells must be ex- 
tracted. The location of the intersections on the cells must 
be found in order to calculate Ssi in Eq. || As mentioned 
before, the searching for the intersecting cells and the val- 
ues stored in the cells are relatively fast since the data is 
stored in the tree data structure. 

Finding the intersections of a photon beam with model 
grids is relatively simple for a spherically symmetric sys- 
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Fig. 6. This diagram shows a photon at r moving in di- 
rection N within the model space. The beam of the photon 
intersects with different sizes of cells before it reaches the 
outer boundary. The value of the opacity of each cell is 
assigned at the center of the cell, and the optical depth 
is calculated simply as the sum of the product of the cell 
value (Xi) and the line segment length (5si) as in Eq. |J. 

tern. In the tree grid case, the following must be done to 
find the intersecting points: 

1. For all the planes of the cells which could have inter- 
sected with this photon, the following set of equations 
must be solved for (x,y,z). The first one is the equa- 
tion of the line along which the photon would move, 
and the second is the equation of the plane. 

(g - -go) = (v - yo) _ ( z - go) /gs 

N x N y N z [ ' 

(x - x c ) n x + (y- y c ) n y + (z - z c ) n z = (6) 

where {x c ,y c ,z c ) is a point on a face of the cube, 
(xo, yo, Zq) is the current location of photon, n = 
(n x ,n y ,n z ) is the normal vector of a face, and N = 
(N x , N y , N z ) is the directional vector of the photon. 

2. There may be up to three such intersections (since the 
ranges of the coordinates are not restricted in Eqs. || 
and |^), but we are only interested in the first inter- 
section which a photon encounters. The distance to 
each intersection must be calculated, and the smallest 
distance is used to compute Ssi in Eq. 

Fig. |^ shows a simple example of how to calculate the 
optical depth between the current photon location (r) and 
the outer boundary (the largest cell in the figure). The 
figure is again depicted in 2-D for clarity, but the same 
method can be used for the 3-D case. The thick solid line 
indicates the cell of current interest. The corresponding 
steps in the figure are the following: 
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1. There is only one intersection of the ray and the cur- 
rent cell (thick line). Since this is not a leaf cell (the 
lowest level cell), descend to the children cells. 

2. There is one intersection, and this is a leaf cell. Find 
the distance (5si) between the photon position (r) and 
the intersection. Compute dr\ = Xl$ s i- Save dr±. 

3. There are two intersections, and this is a leaf cell. 
Find the distance (5s 2) between the two intersections. 
Compute g?T2 = X2^S2. Save dr-i- 

4. There is no intersection, and this is not a leaf cell. 
Dismiss all the subcells. 

5. There are two intersections, but this is not a leaf cell. 
Descend to the children nodes. 

6. There are two intersections, and this is a leaf cell. 
Find the distance (5 S3) between the two intersections. 
Compute drs = X3<5s3. Save dr^. 

7. There are two intersections, and this is a leaf cell. 
Find the distance (5s 4) between the two intersections. 
Compute dT4 — X4<5 S 4- Save c?T4. 

8. There is no intersection, and this is a leaf cell. Do noth- 
ing. 

9. There are two intersections, and this is a leaf cell. 
Find the distance between (5s 5) the two intersections. 
Compute drs = X5^ s 5- Save dT$. 

After the operation described above is finished, the 
summation in Eq. || is performed to find the optical depth 
between the current position of the photon and the outer 
boundary. 

If the cell sizes and their locations are chosen prop- 
erly, the accuracy of the optical depth value is within a 
few percent in normal runs of our models. This method 
is used because of the simplicity and the speed. On the 
other hand, if accuracy of the optical depth value is cru- 
cial, we can use neighboring cell values to improve the 
optical depth estimate. Computing the optical depth is 
the most time consuming part of the code. 

3. Tests 

3.1. Continuum Polarization Tests 

Firstly, the continuum polarization arising from an opti- 
cally thin envelope of a single star is examined. The light 
source is assumed to be a point source which is embedded 
in the density of the following form. 



p(r, 9) = Po-^a 



(1 + 6 cos 2 9) 



(7) 



where p is a constant, a is an angular normalization con- 
stant and b is the parameter which controls the angular 
dependency of p. When b is < 0, = 0, and > 0, the stel- 
lar envelope is oblate, spherical, and prolate respectively. 
The opacity of the atmosphere is set to be r « 0.1. In 
Fig. ||, the normalized polarization levels computed as a 
function of inclination angle are shown, and the y are com- 
pared with the predictions from the model of Brown & 
McLean (1977| ). b = -0.5 is used in both models. The 
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Fig. 7. An illustration of how the optical depth between 
the outer boundary (the largest cell) and the current lo- 
cation (r) of the photon moving in the direction N is cal- 
culated. The cor resp onding instructions for each step are 
summarized in § |2.7| . The thick solid lines indicate the im- 
portant cell in each step. The filled circle is the current 
location the photon, and the open circles are the intersec- 
tions of the photon beam and the cell of current interest. 



figure shows that the results from the two models are in 
excellent agreement. 

Secondly, a binary system which consists of two identi- 
cal point sources orbiting in a circular orbit is considered. 
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Inclination (i) in degrees phase 



Fig. 8. Comparison of the polarization v.s. inclination 
angle of a oblate atmosphere around a single star accord- 



Fig. 9. Comparison of Q & U polarization double sine 



ing to the model of Brown fc McLean (1977j) (lines) and 



our model (circles). The model consist of a point source 
embedded in the density in the form of Eq. [?]. The result 
from our model agrees with the analytic solution of Brown 
& Mcl jean (1977| ) very well. 



Around each light source, a uniform spherical density is 
placed. Again, the envelopes are set to be optically thin 
so that the variation of polarization level can be com- 
pared with the binary polarization model of Brown et al 
(1978[j] Fig. H shows that the results from the two mod- 
els agree each other very well for three inclination angles 
(i = 0°, 45°, 90°). 

Thirdly, we verified the computati ons against t he 2- 
D Monte Carlo polarization model of Hillier (1991 ) and 
the 2-D radiative transfer codes of Hillier (1994 ) whose 
models allow for multiple scatterings (optically thick) and 
an extended source. Our results are in good agreement 
with both models i.e., / and Q are within a fractional 
error of 0.01. 



3.2. Line Polarization Tests 

A He I (A5411) line emitted from a single WN type star 
(R = 1.5%), L = 2.0 x 10 5 Lq, M = 0.6 x 10" 5 M yr" 1 ) 
is computed with the 3-D M onte Carlo co de and the 2- 
D radiative transfer codes of Hillier (1996 ). The angular 



dependency of the atmosphere is assumed to be the same 
as in Eq. [?] with b = —0.5. Fig. [l(] shows the results from 
both models are in excellent agreement. 

An example of a non-axisymmetric system is a rotat- 
ing atmosphere of a single star. Although density distri- 
bution is usually symmetric around a rotation axis, the 
velocity field is not axisymmetric. To examine a basic be- 
havior of polarization variations across an emission line, 
the following two simple models are computed. 1. an ex- 
panding atmosphere (beta velocity law with (3 = 0.5, 



curves predicted by the analytic model of Brown et al 
(1978 ) and our models for three different inclination 



angles. Two identical point sources, which are orbiting 
around each other in a circular orbit, are each embed- 
ded in an optically thin spherical envelope. The points 
with makers are from our model, and the lines are from 
the model of Brown et al. (1978): Circles follow i = 0° 



(face-on) curve, squares follow i = 45°, and diamonds fol- 
low i = 90° (edge-on) curve. The agreement of the two 
models is excellent for all inclination angles. The small 
deviations of our model from that of Brown et al. (1978| ) 
around phase = 0.42 and 0.58 for i = 90° is due to the at- 



mospheric eclipse since the model of Brown et al. (1978| ) 



did not include the attenuation of the star light before 
scatterings. 



Voo — 500 km s" 1 ) and 2. an expanding & rotating atmo- 
sphere (beta- velocity law+solid-body rotation with V m t = 
182 km s -1 on the stellar surface). The results are shown 
in Fig. [ll]. Qualitatively, the behavior of these line varia- 
tions are the same as those described by McLean (1979 ) 
(in their Figure 5), an d those of a simple analytical model 
by |Wood et al. (1993|) (in their Figs. 4, 6 and 8). 



4. Example Calculations for a Binary System: 
V444 Cyg 

As an application to a real system, the phase dependent 
Stokes parameters (/, Q, U) of the massive binary sys- 
tem V444 Cyg (WN5+06 III-V) are computed with our 
model. The system is a short-period (P — 4.212 days, 



Khaliullin ct al. 1984) eclipsing binary, and it exhibits 



variability in polarization, line strength and X-ray flux as 
a function of orbital phase. The variability arises from oc- 
cultation of the photosphere, from perturbations induced 
in the extended atmosphere of the W-R star by the O star 
and its wind, and from the wind-wind interaction region. 
Despite the complexities, many authors (e.g., Hamann 
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Fig. 10. A model of an optically thick Hen (A5411) 
line for a WN5 type star. Normalized flux (top), percent- 
age polarization (middle) and polarization angle (bottom) 
were computed by our 3-D Monte Carlo model (solid) , and 
compared with the 2-D radiative transfer model (dotted) 
of Hillier (1996 ). The angular distribution of electron gas 
is assumed to be the same as in Eq. @ with b = -0.5. The 
agreement between the two models is excellent. 



& Schwarz 1992; St-Louis et al. 1993; Marchenko et al 



; |phcrcpashchuk et al. 1995 




Moffat & Marchenko 


; l^larchcnko et al. 1997 


; |Stevens & Howarth 1999, 



and etc.) have used this object to determine fundamental 
parameters of the W-R star by taking advantage of its 
variable nature. 

The model consists of three main components: 1. the 
W-R star atmosphere, 2. the O star's spherical sur- 
face (with a limb-darkening law) and 3. the paraboloid- 
shaped bow shock region due to the colliding stellar winds 
(Fig. Ilj). More detailed discussion on the V444 Cyg polar- 



ization model can be found in Kurosawa fc Hillier (2001 ) . 
Fig. [l3] shows the resulting continuum Stokes parameters 
(/, Q, U) at 5630A as a function of orbital phase. The 
model reproduces the eclipsing / light curve, the dou- 
ble sine curves of Q & U polarizations and the small 
oscillations seen in Q & U near the secondary eclipse 
(phase=0.5). Next, the phase dependent Hei (A5876) 
emission line, created mainly in the outer part of the W- 
R atmosphere in the binary, is computed. The results are 
shown in Fig. The left plot in Fig. [l4| shows the se- 
quence of He 1 (A5876) line profiles computed with un- 
realistically high emissivity in the bow shock region, to 
emphasize the effect of the bow shock. The underlying or- 
dinary He 1 emission from the W-R star is barely visible 
in this plot. On the left in Fig. [li], the sequence of the 
same line is displayed, but with a more realistic amount 
of emission from the bow shock region. This figure qual- 
itatively describes th e variability seen in th e observation 
of this line (see e.g., Marchenko et al. 1994 ). 




0.00 
-0.20 

0.015 

0.000 



-4000 -2000 2000 4000 

V [km/sl 

Fig. 11. The effect of rotation on line polarization is 
demonstrated in this diagram. Solid line: an emission line 
produced by a radially expanding atmosphere. Dotted 
line: same emission line, but solid-body rotation is added 
to the expanding atmosphere. An inclination angle of 60° 
was used for both models. The angular distribution of 
electrons was assumed to be the same as in Eq. ^ with 
b = -0.5. 




Fig. 12. This figure illustrates the model configuration 
of V444 Cyg (WN5+06 III-V) on the orbital plane. The 
W-R star is placed at the center of the cubic boundary, 
and it is surrounded by the spherical atmosphere. The O 
star is located below the W-R star in the diagram, and the 
tilted paraboloid shock with a given thickness is covering 
the O star. The strong stellar wind from the W-R star is 
interrupted by the shock front, and the density behind the 
shock is assumed to be insignificantly small for simplicity. 
The direction to an observer at phase = 0, 0.25, 0.5, 0.75 
are indicated at the top, bottom, left and right edges since 
the orbital inclination is about 80° (almost edge-on view). 
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0.75 

phase 



Fig. 13. The figure shows the model calculation of polar- 
ization light curves for V444 Cyg. The relative flux (top), 
Q (middle) and U (bottom) polarization at A = 5630A, 
as a function of binary phase, are plotted. Circle: obser- 
vation, Solid: our model with M = 0.6 x 10 -5 M^yr^ 1 . 
The optical light curve data and optical polarization data 
are from Kron fc Gordon (1943 ) and [gt-Louis ct al. (1993 ) 
respectively. 





5. Conclusions 

We have presented a new 3-D Monte Carlo model which 
computes the variability in line and continuum polariza- 
tion associated with the orbital motion of a binary sys- 
tem surrounded by a non-axisymmetric envelope. The ba- 
sic model is constructed in a general manner so that it 
can handle an arbitrary distribution of gas density. A 
special method of grid spacing called the "cell-splitting" 
method is used to automatically assign more grid points 
to the higher density regions. For a fixed number of grid 
points, the tree method achieves a higher accuracy than 
the uniform cubic grid method. This kind of grid scheme 
is very useful in multi-dimensional calculations. It is not 
restricted to a binary system, but also applicable to many 
other astrophysical systems, e.g., accretion flows, molecu- 
lar clouds and a 3-D hydrodynamical model. 

The model was tested extensively. We demonstrated 
the accuracy of our model by comparing it with simple 
analytical models and existing well-tested 2-D numerical 
models. Firstly, the model was tested against the analytic 
models of continuum polarization by Brown & McLean 
(1977) and Brown et al. (1978), and they were in cx- 



-3000 -2000 -1000 1000 2000 3000 
V [km/s] 



-3000 -2000 -1000 1000 2000 3000 
V [km/s] 



Fig. 14. Left: Shows the phase dependent emission from 
the bow shock region (see Fig. [f2]) seen on the top of He I 
(A5876) line. The emissivity of the bow shock region is 
set to be unrealistically strong in order to demonstrate 
the effect of the bow shock emission. Right: Same as for 
the left figure, but the emissivity of the bow shock region 
is reduced by a factor of ~ 10. This figure qualitatively 
describes the variability seen in the observation of this 
line (see e.g., Marchenko et al. 1994). The profiles are 



plotted from phase = (bottom) to phase = 1 (top) with 
0.1 phase steps. 



cellent agreement for both a point source and a binary 
model. Secondly, the polarization in the optically thick 
Hen (A5411) line calculated by our model and by the 
2-D radiative transfer model of [Hillier (1996| ) were com- 



pared. They also showed good agreement. Thirdly, the ef- 
fect of rotation on the line polarization seen in an emission 
line was demonstrated. The basic behavior of polarization 
and polarization angle across an emission line, computed 
by our model, are confirmed to be the same as those in 



McLean (1979) and Wood et al. (1993) 



We also demonstrated the application of our model to 
a real system: V444 Cyg. The model was able to produce 
a set of complicated /, Q and U light curves which fit 
the observational data very well. In addition, we quali- 
tatively modeled the phase dependent excessive emission 
seen on the top of the optically thin Hei (A5976) line. 
Since this excessive emission is originated from the bow 
shock heated region due to the colliding stellar winds, by 
modeling the sequence of this line properly, we would be 
able to probe little known properties (e.g., gas flow, ge- 
ometrical configuration) of the bow shock in the system. 
If one wishes to model a line variability associated with 
the bow-shock region correctly, one needs to have realistic 
opacity information in the region. A more formal approach 
to this kind of problem is to develop a full 3-D non-LTE 
radiative transfer model. Because of the high dimensional- 
ity, computational time would increase much more than in 
the case of a 1-D problem. Developing a 3-D Monte Carlo 
non-LTE radiative transfer model would be simpler than 
developing the 3-D short-characteristic ray method e.g., 
by Folini fe Walder (1999 ), and the tree data structure is 
relevant to this problem. 



In a following paper ( Kurosawa fc Hillier 2001 ), we will 



apply our model, in conjunctio n with the multi-line n on- 
LTE radiative transfer model of |ffillicr fc Miller (1998] ), to 



R. Kurosawa and D. J. Hillier: Tree- Structured Grid Polarization Model 



11 



estimate the mass-loss rate of V444 Cyg. This will be done 
by fitting the observed Hci (A5876) and Hen (A5412) line 
profiles, and the continuum light curves of three Stokes 
parameters (I, Q, U) in V band simultaneously. 
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